Periodicity of mass extinctions without an extraterrestrial cause 
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We study a lattice model of a multi-species prey-predator system. Numerical results show that 
for a small mutation rate the model develops irregular long-period oscillatory behavior with sizeable 
changes in a number of species. The periodicity of extinctions on Earth was suggested by Raup and 
Sepkoski but so far is lacking a satisfactory explanation. Our model indicates that this might be a 
natural consequence of the ecosystem dynamics, not the result of any extraterrestrial cause. 



The Earth ecosystem is certainly a subject of inten- 
sive multidisciplinary research. Researchers in this field 
believe that at least some basic understanding of this im- 
mensely complex system can be obtained using relatively 
simple models, that nevertheless grasp some aspects of 
its rich behavior Q. Of particular interest in physicists 
community is the dynamics of extinctions of species 
Palaeontological data, that show broad distributions of 
these events in the Earth history, suggest existence of 
strong, perhaps power-law correlations between extinc- 
tions. Similar correlations appear in the so-called critical 
systems and such an analogy resulted in a wealth of inter- 
esting models that consider extinctions as a natural con- 
sequence of the dynamics of an ecosystem However, 
fossil data are not entirely convincing, and it is not clear 
to what extent the analogy with critical systems hold. 
What is more, a number of researchers prefer an alterna- 
tive explanation where extinctions appear due to exter- 
nal stresses imposed on the ecosystem as, e.g., impacts 
of comets or meteorites, or an increased volcanic activ- 
ity . The popularity of theories of exogenous origin of 
extinctions increased when Raup and Sepkoski concluded 
from analyzing fossil date that big extinction events dur- 
ing the last 250 My (million years) have been occurring 
with periodicity of about 26 My 0. Several theories, 
mostly of astronomical origin, have been proposed to ex- 
plain such a periodicity, but none of them is confirmed 
or commonly accepted |3- Although the Raup and Sep- 
koski analysis was put into question [gI, the more recent 
analysis confirms a similar periodicity of extinctions 
keeping this fascinating hypothesis still open. 

Lacking a firm evidence of any exogenous cause, one 
can ask whether the periodicity of extinctions can be ex- 
plained without referring to such a factor. Or in other 
words, if it is possible that the ecosystem dynamics pro- 
duces (by itself) oscillations on such a long time scale. 
Since the seminal work of Lotka and Volterra, an os- 
cillatory behavior is already well-known in various prey- 
predator systems P, 0| , but the periodicity of oscillations 
of densities in such systems, that is determined by the 
growth and death rate coefficients of interacting species, 
is of the order of a few years rather than millions. Prey- 
predator systems, where such an oscillatory behavior was 
studied, are typically quite simple and consist of a fixed 
and rather small number of species. Certainly a model 



capable of describing the dynamics of extinctions should 
include a large number of species as well as mutation 
and competition mechanisms. There is already a wealth 
of papers where various models of this kind where ex- 
amined 9] , but none of them has been reporting a long- 
term periodicity of extinctions. There is, however, one 
aspect that these models are missing and that is per- 
haps quite important, namely they neglect spatial cor- 
relations between organisms. From statistical mechanics 
we already know that when the spatial dimension of the 
embedding space is rather low, such correlations might 
play an important role, and hence more realistic models 
of the ecosystem should take them into account. 

In the present paper we study a multi-species lattice 
model of an ecosystem. In our model predator species 
compete for food (prey) and space (to place an offspring) . 
This competition combined with a mutation mechanism 
leads to the periodic behavior, although some character- 
istics of our model, as, e.g., the number of species, show 
in addition strong stochastic irregularities. Sometimes 
our system is populated by a group of medium-efficiency 
species. But this coexistence at a certain moment is in- 
terrupted by creation of a species that is more efficient 
and able to invade even a substantial part of the system. 
However, the reign of such an apex predator does not 
last long. It is a fast-consuming species and it quickly 
decimates the population of preys, which in turn leads 
to its own decline. Such a situation opens up niches that 
again become occupied by less-effective species that sur- 
vived the invasion or were created by mutation, and the 
situation repeats. Simulations show that the smaller the 
mutation probability, the larger the periodicity of such 
a behavior. Although it is difficult to access, we expect 
that the mutation rate in real ecosystems, as interpreted 
in the context of our model is very small and the 
presented model might at least suggest an explanation of 
the 26My periodicity of big extinctions as a natural con- 
sequence of the ecosystem dynamics, not as the result of 
an external perturbation. 

Our model is a multi-species extension of an already 
examined prey-predator model 0. At each site i of 
a square lattice of linear size N we have the four-state 
operator Xi that corresponds to this site being empty 
{xi — 0), occupied by a prey {xi — 1), by a predator 
{xi = 2), or by both of them {xi = 3). Each preda- 
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tor is characterized by a real number parameter rui 
(0 < rui < 1) that we will call size (m^ is meaningful only 
when i is occupied by a predator). We also introduce the 
relative update rate of preys and predators r (0 < r < 1), 
and the mutation probability p. The dynamics of this 
model is specified as follows: 

(a) Choose a site at random (the chosen site will be 
denoted hy i). 

(b) With the probability r update a prey at site i (i.e., 
if Xi = 1 or Xi = 3, otherwise do nothing). Provided 
that at least one neighbor (say j) of the chosen site is 
not occupied by a prey (i.e., Xj = or Xj = 2), the prey 
at the site i produces an offspring and places it on an 
empty neighboring site (if there are more empty sites, 
one of them is chosen randomly). Otherwise (i.e., if 
there are no empty sites) the prey does not breed. 

(c) Provided that i is occupied by a predator (i.e., 
xi — 2 or Xi = 3) update this site with the probability 
(1 — r)mi, where rrii is the size of the predator at site i. 
If a chosen site is occupied by a predator only (xi = 2), 
it dies, i.e., the site becomes empty {xi = 0). If there 
is also a prey there {xi = 3) the predator consumes the 
prey (i.e., Xi is set to 2) and if possible it places an 
offspring at an empty neighboring site. For a predator 
of size rrii it is possible to place an offspring at a site j 
provided that j is not occupied by a predator [xj = 
or Xj = 1) or is occupied by a predator [xj = 2 or 
Xj = 3) but of a smaller size than rrii (in such a case 
a smaller-size predator is replaced by an offspring of a 
larger-size predator). An offspring inherits parent's size 
with the probability 1 — p and with the probability p it 
gets a new size that is drawn from a uniform distribution. 

One can see that the size rrii of a predator deter- 
mines both its update rate and its strength when it 
competes with other predators for space. While the in- 
creased strength is always favorable, the larger update 
rate might be a disadvantage when preys do not repro- 
duce fast enough. As it will be shown below, the behavior 
of our model is very much influenced by this property of 
the dynamics. 

The already studied single-predator version is ob- 
tained when all predators have a unit size rui = 1 and 
suppressed mutations p = 0. In such a case, for r > 0.11 
the model is in an active phase with positive densities 
of preys po (which is a fraction of all sites i such that 

= 1 or Xi — 3) and predators p (fraction of all sites 
i such that Xi = 2 or Xi — 3). For r < 0.11 the update 
rate of preys is too small to sustain an active phase but 
it is a population of predators that becomes extinct and 
the model enters an absorbing state where all sites are 
occupied by preys. In the active phase but close to the 
transition point (0.11) one observes oscillations of po and 
p but the amplitude of these oscillations diminishes in 
the thermodynamic limit N ^ oo. On the other hand, 
for the model on the three-dimensional lattice such oscil- 
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FIG. 1; Average size m, fraction of a dominant predator 
species / and the number of species s as a function of update 
rate r. Results of simulations do not depend on an initial con- 
figuration and usually it was a random distribution of preys 
and predators. 



lations most likely persist in this limit llll |. 

To examine the behavior of our model we used simu- 
lations and measured its various characteristics such as 
densities of preys po and predators p, the average size of 
dominant predator /, the average size to, the number of 
species s, and the lifetime of predator species. To define s 
we classify predators into species according to their size. 
Some of these quantities are presented in Fig. One 
can see that for r > Vc ^ 0.27 predators in the system 
belong essentially to one dominant (/ ~ 1) species of a 
large size (m ^1). Of course, mutations create from 
time to time some other species but they occupy a negli- 
gible portion of a system - unless a newly created species 
will have a larger size than the dominant species and will 
be able to invade the system. Fig. ^ also shows that a 
much different behavior appears for r < re- In this case 
a dominant species occupies only a small fraction of a 
system (the comparison with the results for system size 
N — 200 shows a strong iV-dependence and suggests that 
for larger N the fraction / will diminish to zero). More- 
over, the average size m differs substantially from unity 
that indicates that having a large size is no longer advan- 
tageous feature. Another indication of a more complex 
behavior in this case is a large increase of the number of 
species. 

In our opinion, it is the regime for r < Tj, whose com- 
plex dynamics might resemble the behavior of realistic 
ecosystems. To have a better understanding of the be- 
havior of the model in this regime we present a time de- 
pendence of some of its characteristics. The unit of time 
is defined as a single, on average, update of each site (i.e., 
it is made of N'^ elementary single-site updates). While 
in Fig. [3 densities po, p, and the average size to show 
a relatively regular oscillations, the number of species s 
is much more irregular. During periods of multi-species 
coexistence, predators have a rather small size (they eat 
slowly) that enables them to sustain their density p rel- 
atively large. As a result density of preys po is rather 
small. At certain moment, however, a predator of a large 
size is created and starts to invade the system. As a 
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FIG. 2: Time dependence of the number of species s (to su- 
perpose with other data it was divided by 40), average size 
m, density of preys po, and density of predators p. 



result the number of species s rapidly decreases while 
771 increases. Moreover, the density p decreases and this 
is related with the fact that a predator of a large size 
consumes preys too quickly and is simply running out 
of food. Hence, the population of this predator in some 
places disappears and that creates areas where preys can 
breed without being consumed by predators and that is 
why the density po after an initial short decline increases 
to a relatively large value. However, a large-size species 
cannot keep its dominance for a long time since large 
empty places occupied mainly by preys constitute ideal 
niches for other predators as well. As a result, the model 
is driven again toward a multispecies coexistence. 

An important question is how these oscillations be- 
have for an increasing system size N. Comparing (not 
presented) results for different values of N, we expect 
that the amplitude of these oscillations will diminish to 
zero (period of oscillations does not seem to depend on 
N) . This is because for a sufficiently large N the system 
is essentially decomposed into several independent do- 
mains where multi-species and fewer-species periods are 
uncorrelated and fluctuations cancel out. However, there 
is an additional factor that is responsible for the size of 
these independent domains and thus the amplitude of 
oscillations, namely the mutation probability p. Indeed, 
the end of the multi-species period in a certain domain 
is induced by the creation of a large-size predator. For 
the decreasing mutation probability p such events will 
be less and less frequent and multi-species domains will 
have more time to grow. We thus expect that for decreas- 
ing p the size of such domains should increase and, as a 
result, for finite N the amplitude of oscillations should 
also increase. Moreover, the period of these oscillations, 
that is determined by the time needed for such domains 
to grow, should also increase. Simulations, as shown in 



FIG. 3: The time dependence of the number of species for 
(from top) p = 0.00001, 0.0001, 0.001, and 0.01. To super- 
pose the data on a graph the actual values of s were divided 
by some factors. Such an operation does not change a char- 
acteristic period of fluctuations and their relative amplitude. 
Inset shows the period of oscillation r as a function of muta- 
tion probability p obtained from the maximum of the Fourier 
transform of the time dependence of the number of species 
{N = 1000) 

Fig. 13 confirm such a behavior. Let us notice large fluc- 
tuations for p ~ 0.00001, where the number of species 
after an invasion drops roughly by a factor of two. To 
examine the p-dependence of the period of oscillations t 
more quantitatively, we calculated the Fourier transform 
of the time dependent number of species s (other charac- 
teristics like 777, pq, or p give basically the same result). 
The period of oscillations t as extracted from the maxi- 
mum of this transform is shown on the logarithmic scale 
in the inset of Fig. 13 Straight line fit corresponds to the 
dependence t p-o.3i calculations for larger system 
size N or smaller p might modify this estimation. 

As we already mentioned, the amplitude of oscillations 
in our model is determined by the combination of two 
factors: system size N and mutation probability p. That 
T increases for decreasing p is an important result. It 
shows that for a small mutation probability p and fi- 
nite but large system size N (i.e., specifications of the 
real ecosystem) the model develops long-period oscilla- 
tory behavior with sizeable changes of e.g., the number 
of species s. It might be interesting to notice that for 
the single-predator version |lll |. with rrii — 1 and p = 0, 
the period of oscillations in the two-dimensional case for 
r = 0.2 is around 30. For the present model and for 
p = 0.00001 the period of oscillations is larger by almost 
three orders of magnitude (see the inset in Fig. OJ. It 
shows that the periodic behavior in our model has a much 
different mechanism than the Lotka-Volterra oscillations 
in simple prey-predator systems. 

One of the properties often analyzed in models 
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FIG. 4: Logarithmic plot of the probabihty distribution of 
hfetimes of predator species. 

of ecosystems is the lifetime distribution of species. 
Palaeontological data suggest some broad distributions 
but they are again not very conclusive and both expo- 
nential and power-law fits can be made . The lifetime 
distribution for our model is shown in Fig. ^ Although 
for p = 0.01 the numerical results suggest an exponen- 
tial distribution, for smaller p the situation is less clear. 
Especially, for p = 0.00001 it seems that a broader, per- 
haps a power-law distribution might better describe the 
lifetime of our species. 

It would be interesting to make further analysis of our 
model. For example, one can implement a less abrupt 
mutation mechanism, where a new species will be only 
a small mutation of a parental species. Such a modifi- 
cation probably results in a longer period of oscillations 
and might be more suitable for comparison with the real 
ecosystem. Another possibility might be to examine the 
differences in, e.g., lifetime distribution of species before 
and after an extinction. That such differences exist is 
suggested by the asymmetry of our data in Fig. |51 where 
the changes in a pre-extinction period seem to be dif- 
ferent than in the post-extinction one. Palaeontological 
data also show certain differences in longevity of species 
during such periods , and a comparison with the pre- 
dictions of our model, if feasible, would be very desirable. 
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